#ifndef ATMOSLAYER_H
#define ATMOSLAYER_H

#include <QObject>
#include <QMap>
#include <QByteArray>
#include "mie_c_src/complex.h"
#include <QList>
#include <math.h>

//Stokes parameter
#define I 0
#define Q 1
#define U 2

typedef QMap<long int, double> H_NUM;
typedef QMap<double, double> WL_CROSS;


#define N_PI 90 //zenith and azimuth angle partition numbers for Pi
#define N_ZENITH N_PI/2
#define ANGLE_STEP (M_PI/N_PI)
#define ZENITH_INIT (ANGLE_STEP/2)

#define N_SCAT 3001 //scatter angle partition
#define N_SMALL_SCAT 2001
#define SMALL_SCAT_REGION 0.1
#define STEP_SMALL_SCAT (SMALL_SCAT_REGION/(N_SMALL_SCAT-1))
#define STEP_BIG_SCAT ((M_PI-SMALL_SCAT_REGION)/(N_SCAT - N_SMALL_SCAT))

#define DOWN true
#define UP false

#define UP_IN_UP_SCAT 1
#define UP_IN_D_SCAT 2
#define D_IN_UP_SCAT 3
#define D_IN_D_SCAT 4

#define SAVE_DISABLE 0x00
#define SAVE_UP 0x01
#define SAVE_DOWN 0x10
#define SAVE_ALL 0x11


#define Sin sin
#define Cos cos
#define Power pow
#define Sqrt sqrt

typedef struct {
	double i, q, u;
}linearPolar ;

/*! \~english
    \class AtmosLayer
    \brief Class AtmosLayer is the base class in VRTCaculate library.

    This class is cross-platform using Qt4 for peripheral(file read/write) and  container classes(list, vector...), and can be used to create applications Qt4.

    The most important role class AtmosLayer plays is that it defines the interface for any kind of user realized atmospheres(Rayleigh, Aerosol..), which is:
    <TABLE>
    <TR><TH> Name </TH><TH> Signal/Slots </TH><TH> Note </TH></TR>
    <TR><TD> emitSolarLight </TD><TD> signal </TD><TD> Emit solar light to lower layer, stimulate recvSolar. </TD></TR>
    <TR><TD> recvSolar </TD><TD> slot </TD><TD>  Recieve solar irradiance and direction.</TD></TR>

    <TR><TD> emitLightUp </TD><TD> signal </TD><TD> Emit upward scattering light to upper layer, stimulate recvScatInUp. </TD></TR>
    <TR><TD> recvScatInUp </TD><TD> slot </TD><TD>  Recieve and cache upward incident scattering light(from bottom of layer), zenith (Pi/2, Pi)</TD></TR>

    <TR><TD> emitLightDown </TD><TD> signal </TD><TD> Emit downward scattering light to lower layer, stimulate recvScatInDown.  </TD></TR>
    <TR><TD> recvScatInDown </TD><TD> slot  </TD><TD>  Recieve and cache downward incident scattering light(from top of layer), zenith (0, Pi/2)</TD></TR>
    </TABLE>

    For middle placed(default) layers, the flow path is:\n
    \b Upward:
    emitLightUp(lower) -> \n
    recvScatInUp(upper, cache && start a 10ms timer && return) ->\n
    upInResponse -> upInUpScat() and upInDownScat().

    \b Downward: \n
    emitSolarLight(upper, if has solar light) -> \n
    recvSolar(lower, cache && return) -> \n
    emitLightDown(upper) -> \n
    recvScatInDown(lower, cache && start a 10ms timer && return) ->\n
    downInResponse ->\n
    solarScatCaculate -> solarUpScatCaculate() and solarDownScatCaculate() ->\n
    downScatInResponse -> downInDownScat() and downInUpScat()\n
\note For upper layers, signal emitSolarLight must be emit first(if has solar light), then emitLightDown, because slot recvScatInDown only cache solar parameter, recvScatInDown call all the scattering-caculate process.

    These interface function can be virtual-inherit to achieve specific behavior of top/ground layer.\n
    \b TopLayer:
    There is no scattering incident light at top of earth atmosphere, so there can be no emitLightDown signal binded at recvScatInDown for a TopLayer. So recvSolar() must be modified(virtual-inherit) to call downInResponse(), who is changed accordingly(no scatter incident response caculate). And, recvScatInUp() must be modified, that is, when caculating upward incident light's responce, it must change caculating direction, and emitLightDown.

    \b GroundLayer:
    Similarly, it change direction when recieved downward light, and emit ground scattering light according to BRDF, and finally emit groundScatUp, which get recvScatInUp() called.

    Precaculated scattering-template is used to speed up caculation. Scattering-template is stored in template repository, one resository is compose of series of .tmpl files for various optical depth(tau) for the same kind of atmosphere(phase function), who is specified by a file "atmos.spec.xml", which can be searched for later use.
 */

class AtmosLayer : public QObject
{
	Q_OBJECT
public:
	explicit AtmosLayer(QString layerName,
			    double hBottom,
			    double hTop,
			    QObject *parent = 0);
	~AtmosLayer();
/*****************Interface: out************************************/
signals:
    void emitLightUp(QByteArray lightData);
    void emitLightDown(QByteArray lightData);
    void emitSolarLight(double intensity,
			double zenith, double azimuth);
    void gotDownIrradiance(double irradiance);

protected:
    void emitLight();

/****************Interface: in(see html documents for call sequence)***/
public slots:
    virtual void recvSolar(double intensity,
			   double zenith, double azimuth);
    virtual void recvScatInUp(QByteArray lightData);
    void recvScatInDown(QByteArray lightData);


/***************Kernel: caculating stuff****************************/
protected slots: //called by a asynchronous timer to allowed SIGNAL(emitLight) return.
    virtual void upInResponse();
    virtual void downInResponse();

protected:
    //downward incident has 2 parts, scatter && solar, get called in downInResponse().
    virtual void solarScatCaculate();
    virtual void downScatInResponse();

protected: //functions correspond to caculate procedure
    //these 4 is used to speed up caculation, e.g. cache scat-coeff.
    //default behevior is call scatCoeff() every time.
    virtual void upInUpScat();
    virtual void upInDownScat();
    virtual void downInDownScat();
    virtual void downInUpScat();

    //solar scattering caculation is less, no need for optimize generally.
    virtual void solarUpScatCaculate();
    virtual void solarDownScatCaculate();

    //You Gotta Realize these 2 func for a particular kind of phase-func.
    virtual void solarScatCoeff(int iThetaS, int iPhiS,
				linearPolar *lightOut,
				int scatMode);
    virtual void scatCoeff(int iThetaS, int iPhiS,
			   int iThetaI, int iPhiI,
			   double lightIncident[3],
			   linearPolar *lightOut,
			   int scatMode);

protected:
    //parallel incident light parameters: solar intensity, and its directions
    double m_solarIntensity, m_solarZen, m_solarAzi, m_solarOut;

    double (*m_lightInUp)[2*N_PI][3]; //incident upward light at all azimuth
    double (*m_lightInD)[2*N_PI][3]; //incident downward light at all azimuth

    double (*m_lightScatUp)[2*N_PI][3]; //scattered upward light at all azimuth
    double (*m_lightScatD)[2*N_PI][3];//scattered downward light at all azimuth

protected: //******** Result output*****************
    double (*m_lightScatUpTotal)[2*N_PI][3]; //scattered upward light at all azimuth
    double (*m_lightScatDTotal)[2*N_PI][3];//scattered downward light at all azimuth

    void saveResult(bool dir);
    static QString m_resultPath;


/*********************irradiance caculate stuff********************/
protected slots:
    void downIrradianceCaculate();//integrate down intensity moduled by mu
protected:
    double m_downIrradiance;
    double m_solarIrradiance;
    double m_currentIrradiance;


//****************Interface: config ************************
public:
    static void setWaveLength(double wl); //in micron(um, not nm)
    QString getLayerName();

    void setTau(double tau);
    void setResultSaveMode(char mode);
    void setSavePath(QString path);
    QString getResultPath();
    static QList<double> expHeightProfileCaculate(double stdHeight, double lowerH, double upperH, int numLayers);

    static void scatAngleInit();

protected:
    QString m_layerName;
    static bool m_currentDir;
    int m_loopUpTime, m_loopDownTime;

    char m_saveMode;

protected:
    static double m_wl;
    double m_hBottom, m_hTop;

    double m_tau; //optical depth at Z-axis
    double m_moleTau;

    static H_NUM m_moleNumber;
    static H_NUM m_moleDensity;
    void moleTauCaculate(); //caculate molecule optical depth
    double getMoleTauWholeAtmos(); //according to wavelength

protected:
    static double m_solidAngle[N_PI/2]; //solid angle at zenith angle [0, pi/2]
    static void solidAngleCaculate(); //caculate solid angle at N zenith angle

    static QList<double> m_zenithUpList;
    static QList<double> m_zenithDList;
    static QList<double> m_aziList;
    static void zenithAngleInit();

    //iPhiI symmetry: scatAngles(thetaI, thetaS, phiI, phiS) = scatAngle(thetaI, thetaS, 0, phiS - phiI)
    //                        zenIn     zenOut    aziOut
    static double (*m_scatAnglesSameDir)[N_ZENITH][2*N_PI];
    static double (*m_scatAnglesDiffDir)[N_ZENITH][2*N_PI];

protected: //intialization funtions
    void layerInit();

};
//(360×3×3 + 90 + 90×180×3×6)×8÷1024÷1024
#endif // ATMOSLAYER_H
